In this analysis, we use R to read in LandSat8 data, create a vegitation index and calculate zonal statistics of a raster. The data come from the 2006 National Land Cover Database.
First we relcassify the NLCD data into two classes, based on the value of the raster. In this case, I want to classify the pixels as to whether they are developed or undeveloped.
We then use a census tract polygon layer to calculate the proportion of each tract’s land area that is developed vs undeveloped.
library(rgdal)
## Loading required package: sp
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-28, (SVN revision 1158)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Program Files/R/R-4.1.2/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Program Files/R/R-4.1.2/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Program Files/R/R-4.1.2/library/rgdal/proj
library(raster)
library(gdalUtils)
nlcd<-raster("C:/Users/ozd504//OneDrive - University of Texas at San Antonio/classes/dem7093/dem7093_22/data/nlcd16_48_lc/nlcd_2016_land_cover_48_20190424.img.vat.img")
nlcdproj<-gdalwarp(srcfile ="C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/nlcd16_48_lc/nlcd_2016_land_cover_change_index_48_20190424.img.vat.img",
dstfile = "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/outputnlcd.tif",
t_srs = nlcd,
output_Raster = T)
## Warning in system(cmd, intern = TRUE): running command '"C:\Program
## Files\QGIS 3.16.16\bin\gdalwarp.exe" -t_srs "raster(ncols=44253,
## nrows=41260, xmn=-1075365, xmx=252225, ymn=309855, ymx=1547655,
## crs='+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0
## +y_0=0 +datum=WGS84 +units=m +no_defs')" -of "GTiff" "C:/Users/ozd504/
## OneDrive - University of Texas at San Antonio/gis_classwork/nlcd16_48_lc/
## nlcd_2016_land_cover_change_index_48_20190424.img.vat.img" "C:/Users/ozd504/
## OneDrive - University of Texas at San Antonio/gis_classwork/outputnlcd.tif"' had
## status 1
raster::plot(nlcdproj,
main="Original Raster image")
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
##
## Attaching package: 'sf'
## The following object is masked from 'package:gdalUtils':
##
## gdal_rasterize
library(tigris)
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
library(terra)
## terra 1.5.21
##
## Attaching package: 'terra'
## The following object is masked from 'package:rgdal':
##
## project
options(tigris_class = "sf")
counties<-tracts(state="TX",
county = c("Bexar", "Comal"),
cb=T,
year=2016)
counties<-st_transform(counties, crs=2278)
# doesn't work
est<-terra::crop(nlcdproj, counties)
plot(est)
est2<-terra::mask(est, counties)
plot(est2)
# est<-mask(x=nlcdproj,mask=as_Spatial(counties))
# est2<-crop(x=nlcdproj, y=extent(as_Spatial(counties)))
# plot(est2)
m<-c(1, 12, 0,
13, 24, 1,
25, 95, 0)
rclmat<-matrix(m,
ncol=3,
byrow = T)
nlcdreclass<-raster::reclassify(x = est2,
rcl = rclmat)
plot(nlcdreclass,
main="Reclassified Raster image \n 1=developed, 0 = undeveloped",
col= gray(100:0 / 100))
library(tigris)
satract<-tracts(state="48",
county="029",
year=2018)
library(sf)
satract<-st_as_sf(satract)
satract<-st_transform(satract,
crs= 2278)
satract<-as(satract, "Spatial")
summary(satract)
Now we find the mean of the binary reclassified data, to estimate the percent of each tract’s area that is developed.
First we reproject the raster to the same coordinate system as the census data. I get the projection string from ^^^ above.
nlcdreclass<-projectRaster(nlcdreclass,
crs="+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
satract$pct_developed<-as.numeric(extract(nlcdreclass, satract, fun=mean)) #this takes a minute or two
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the crs of the
## Raster
satract<-st_as_sf(satract)
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:terra':
##
## arrow
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
##
## intersect, src, union
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
satract<-satract%>%
filter(complete.cases(pct_developed))
tractlines<-st_boundary(satract)
library(tmap)
tm_shape(satract)+
tm_polygons("pct_developed",style="quantile", n=5,legend.hist=T )+
tm_shape(tractlines)+
tm_lines(col="black")+
tm_format("World",
title="San Antonio Developed Land Area",
legend.outside=TRUE )
I used this method in a paper I wrote in 2013, although in that paper I did the diversity of the land cover. That could be done like:
library(vegan)
## Warning: package 'vegan' was built under R version 4.1.3
## Loading required package: permute
## Warning: package 'permute' was built under R version 4.1.3
## Loading required package: lattice
## This is vegan 2.5-7
shannonVegan <- function(x, ...) {
diversity(table(x), index="shannon")
}
#shanVegOut <- focal(nlcdproj, fun=shannon, pad=T)
satract$var_land<-as.numeric(extract(nlcdproj,
satract,
fun=shannonVegan)) #this takes a minute or two
shannon<-tm_shape(satract)+
tm_polygons("var_land",
style="quantile",
n=5,
legend.hist=T )+
tm_shape(tractlines)+
tm_lines(col="black")+
tm_format("World",
title="San Antonio Shannon Diversity of Land Cover",
legend.outside=TRUE )
shannon
I recommend using the ggsave() function to export images
made by ggplot for use in other documents. For tmap, you
can use tmap_save()
#ggsave(m1, path = "", filename = "sa_pct_devel.png", units = "in", width=8, height = 10)
tmap_save(shannon,
filename = "~/OneDrive - University of Texas at San Antonio/classes/gis_classwork/shannon_tract.png",
width = 8,
height = 10,
units = "in")
And you can insert that image into a doc or presentation.
Good simple description of landsat here
library(raster)
library(rgdal)
library(rgeos)
## rgeos version: 0.5-9, (SVN revision 684)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.4-6
## Polygon checking: TRUE
library(RColorBrewer)
# turn off factors
options(stringsAsFactors = FALSE)
landsatstack<- list.files("C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/",
pattern=glob2rx("*B*.TIF$"),
full.names = T)
landsatstack
## [1] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B1.TIF"
## [2] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B10.TIF"
## [3] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B11.TIF"
## [4] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B2.TIF"
## [5] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B3.TIF"
## [6] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B4.TIF"
## [7] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B5.TIF"
## [8] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B6.TIF"
## [9] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B7.TIF"
## [10] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B8.TIF"
## [11] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B9.TIF"
## [12] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_BQA.TIF"
lsb2<-raster(landsatstack[2])
plot(lsb2,col = gray(0:100 / 100))
rm(lsb2) ; gc()
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 4349763 232.4 8465130 452.1 8465130 452.1
## Vcells 26157396 199.6 101976122 778.1 451409981 3444.0
lsstack<-stack(landsatstack[c(1,4:9)])
lsbrick<-brick(lsstack)
names(lsbrick)
## [1] "LC08_L1TP_027040_20200311_20200325_01_T1_B1"
## [2] "LC08_L1TP_027040_20200311_20200325_01_T1_B2"
## [3] "LC08_L1TP_027040_20200311_20200325_01_T1_B3"
## [4] "LC08_L1TP_027040_20200311_20200325_01_T1_B4"
## [5] "LC08_L1TP_027040_20200311_20200325_01_T1_B5"
## [6] "LC08_L1TP_027040_20200311_20200325_01_T1_B6"
## [7] "LC08_L1TP_027040_20200311_20200325_01_T1_B7"
names(lsbrick) <- gsub(pattern = "LC08_L1TP_027040_20200311_20200325_01_T1_",
replacement = "",
names(lsbrick))
names(lsbrick)
## [1] "B1" "B2" "B3" "B4" "B5" "B6" "B7"
ndvi<-(lsbrick[[5]] - lsbrick[[4]]) / (lsbrick[[5]] + lsbrick[[4]])
# plot(ndvi,
# main = "NDVI around San Antonio",
# axes = FALSE, box = FALSE)
library(tmaptools)
tm_shape(ndvi)+
tm_raster(style="quantile",
n=5,
legend.hist=T,
palette =get_brewer_pal("RdYlGn", n = 5) )+
tm_format("World",
title="NDVI - San Antonio Region",
legend.outside=TRUE )+
tm_shape(satract)+
tm_polygons(alpha = .7)
## stars object downsampled to 991 by 1010 cells. See tm_shape manual (argument raster.downsample)
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
hist(ndvi)
library(mapview)
mapview(nlcdreclass)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 9494355
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 9494355 '
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
writeRaster(x= ndvi,
filename = "C:/Users/ozd504//OneDrive - University of Texas at San Antonio/gis_classwork/ndvi1_22.tif", format="GTiff", datatype="INT2S", overwrite=T)
plotRGB(lsbrick,
r = 6, g = 5, b = 2,
stretch = "lin",
axes = TRUE,
main = "RGB composite image\n showing Agricultural land - Landsat Bands 6, 5, 2")
box(col = "white")
plotRGB(lsbrick,
r = 5, g = 4, b = 3,
stretch = "lin",
axes = TRUE,
main = "Color infrared composite image\n Landsat Bands 5, 4, 3")
box(col = "white")